#ifndef __USRDENS
#define __USRDENS
#include "type.h"
#define densitye 1e24
#define ATOM_NUM  79.0
#define ngold 5.89e28
#include "physicalconstants.h"
#include "common.h"
#define UNITS 1E-6
void calc_density(const __Vect2<double>& position, vector<double>& density)
{
    double x = position.member[0];
    double y = position.member[1] - 15 * UNITS;
    __Vect2<double>  pos1(18 * UNITS, 1 * UNITS), pos2(28 * UNITS, 29 * UNITS);
    __Vect2<double>  pos3(24 * UNITS, 14.7 * UNITS), pos4(24.3 * UNITS, 15.3 * UNITS);
    __Vect2<double> center(30e-6, 20e-6);
    __Vect2<double> radius = position - center;
    double rad = 10e-6;
    double ncritt = 1.1e27;
    
    if(position > pos1 && position < pos2)
    {
        density = {5.0 * ncritt, 5.0 / 6.0 * ncritt, 0, 0, 0};
    }

    if((y * y < 5.0625 * UNITS * (x - 24 * UNITS)) || (y * y < -5.0625 * UNITS * (x - 22 * UNITS)))
    {
        density = {0, 0, 0, 0, 0};
    }

    if((position > pos3 && position < pos4))
    {
        density = {20 * ncritt, 0, 0, 0, 20 * ncritt};
    }
    
    
    /***
    if(position > pos1 && position < pos2)
    {
        //density[0] = densitye;
        density = {0, ngold, 0, 0, 0};
    }
    //cout<<"x, y "<<x<<"\t"<<y<<"\t"<<(x - UNITS>= 50 * y*y)<<endl;
    //getchar();
    if(x - 1 * UNITS >= y * y / UNITS)
    {
        density[1] = ngold;
    }
    
    double y1 = 0.267 * x + 10e-6, y2 = - 0.267 * x + 30e-6;
    if(x > 40e-6)
    {
        return;
    }
    if(((y > y2 - 1.5e-6 && y < y2) ||( y < y1 + 1.5e-6 && y > y1)) && x < 35e-6)
    {
        density[0] = 4.0 * ncritt;
        density[4] = 0.666 * ncritt;
        //density = {4.0 * ncritt, 0, 0, 0, 0.666 * ncritt};
    }
    if(((y > y2 && y < y2 + 1.5e-6) || (y < y1 && y > y1 - 1.5e-6)) && x < 40e-6)
    {
        density[0] = 30.0 * ncritt;
        density[1] = 5.0 * ncritt;
        //density = {30.0 * ncritt, 5.0 * ncritt, 0, 0, 0};
    }
    ***/
    
    
    /***
    if((y > 0.25 * x + 10e-6) && (y < -0.25 * x + 30e-6) && (x < 20e-6))
    {
        density = {0, 0, 0, 0, 0};
    }
    if((y < 0.25 * x + 10e-6) && (y > 0.25 * x + 6e-6) && \
            ( y < - 0.25 * x + 34e-6) && ( y > - 0.25 * x + 30e-6) && \
            (x < 20e-6))
    {
        density[4] = ngold * 100;
    }
    if(Mod(radius) < rad)
    {
        density[1] = 10 * densitye;
    }
    ***/
}

#endif



#ifndef __INITSP
#define __INITSP

void init_sp(vector<Particles>& spptcls, vector<int>& ppc)
{
    __Vect3<double> velocity(1e8, 1e8, 0), zeros(0, 0, 0);
    /// particle species struct: {name , sp_id, mass, charge, energy, velocity}
    Particles electrons = {"Electron", 0, e_mass, -e_charge, 0.0, zeros};
    //Particles ions = {"Au", 1, 63.5 * 1837 * e_mass, 6.0 * e_charge, 0.0, zeros};
    Particles proton = {"Proton", 4, 1837 * e_mass, e_charge, 0.0, zeros};
    Particles position = {"Positron", 2, e_mass, e_charge, 0.0, zeros};
    Particles photon = {"Photon", 3, 0, 0, 0, zeros};
    Particles Carbon = {"Carbon", 1, 12 * 1837 * e_mass, 6 * e_charge, 0.0, zeros};
    spptcls.push_back(electrons);
    spptcls.push_back(Carbon);
    spptcls.push_back(position);
    spptcls.push_back(photon);
    spptcls.push_back(proton);
    ppc = {5, 5, 5, 5, 5};
}



#endif
